#define CHECK(x) if(x /= 0) write(logunit,*) 'hdf5 error at line ',__LINE__,' of file ',__FILE__

subroutine parse_paramfile()
  use log
  use myfiles
  use read_parameters
  implicit none
  
  write(logunit,*)' -reading parameters from: ',trim(adjustl(paramfile))
  call flush(logunit)
  call read_parameterfile(trim(adjustl(paramfile)))

  call get_parameter('BaseDir'       ,basedir       ,required=.true.)
  call get_parameter('BaseFile'      ,basefile      ,required=.true.)
  call get_parameter('OutDir'        ,outdir        ,required=.true.)

#ifdef METAL_LINES
  call get_parameter('MetalEmissionPath'     ,MetalEmissionPath     ,required=.true.)
  call get_parameter('MetalLinesList'        ,MetalLinesList        ,required=.true.)	
  call get_parameter('MetalLinesRedshiftList',MetalLinesRedshiftList,required=.true.)	
#endif

#ifdef APEC
  call get_parameter('XrayPath'      ,XrayPath      ,required=.true.)
#endif

#ifdef GALAXEV
  call get_parameter('BC03File'      ,BC03File      ,required=.true.)
  call get_parameter('BC03FilterList',BC03FilterList    ,required=.true.)
  call get_parameter('SunSpectrum'   ,SunSpectrum   ,required=.true.)
  call get_parameter('VegaSpectrum'  ,VegaSpectrum  ,required=.true.)
#endif

  call finalize_parameters(abort_on_unused=.true.,abort_on_default=.false.)
  write(logunit,*)''
  call flush(logunit)

end subroutine parse_paramfile

! *****************************************************************

subroutine read_header_params_constants(isnap)
  use log
  use myfiles
  use params_and_constants
  use hdf5_wrapper
  implicit none

  integer(hid_t)     :: file_id
  integer(hid_t)     :: error
  integer*4          :: isnap
  character(len=512) :: infile_no_ext,LocalSnapshot
  logical            :: imulti

  write(csnap,'(I3.3)')isnap 
  write(logunit,*)'  -reading snapshot header, parameters and constants: isnap = '//trim(adjustl(csnap))
  call flush(logunit)
  
  ! Look for snapshot as single file; if not there look for the .0 etc extension 
  ! Turn off hdf5 error messages
  call h5eset_auto_f(0,error)  
  write(suffix,'(I3.3)')isnap 
  infile_no_ext = trim(adjustl(basedir))//'/snapshot_'//trim(adjustl(suffix))&
       //'/'//trim(adjustl(basefile))//'_'//trim(adjustl(suffix))
  
  ! Test for presence of file - if not there, add on .0 to the filename
  snapfile = trim(adjustl(infile_no_ext))//'.hdf5'

  call h5fopen_f(snapfile,H5F_ACC_RDONLY_F,file_id,error)
  if(error.lt.0) then
     ! Can't open, try .0.hdf5 instead and try again
     snapfile = trim(adjustl(infile_no_ext))//'.0.hdf5'
     call h5fopen_f(snapfile,H5F_ACC_RDONLY_F,file_id,error)
     if(error < 0) then
        write(logunit,*)trim(snapfile)
        call abort('   -Cannot find snapshot in single- or multi-file format')
     else
        write(logunit,*)'   -opening snapshot in multi-file format'
        call flush(logunit)
        imulti = .TRUE.
        call h5fclose_f(file_id,error)
     endif
  else
     write(logunit,*)'   -opening snapshot in single-file format'
     call flush(logunit)
     imulti = .FALSE.
     call h5fclose_f(file_id,error)
  endif
  
  ! Turn on hdf5 error messages again
  call h5eset_auto_f(1,error)

  ! Read the various constants
  call read_constants()
  call read_header()
  call read_parameters()
  call read_units()
  call read_cgs_conversions()
  
  write(logunit,*)'   -done'
  call flush(logunit)

end subroutine read_header_params_constants

! --------------------------------------------------------------------

subroutine read_snapshot(isnap,get_mass)
  use log
  use myfiles
  use particles
  use hdf5_wrapper
  use params_and_constants
  implicit none
  integer :: file_handle,rank,dims(8)
  integer :: isnap,ifile,ntotal,strlen
  integer :: i
  integer start2D(2), count2D(2), start1D(1), count1D(1)
  integer icount,ioffset(0:5),ipart
  integer(hid_t) :: file_id,error
  real*4, allocatable, dimension(:,:) :: real2d
  character(len=512) :: PartTypeName,varname
  character(len=512) :: infile_no_ext,LocalSnapshot
  logical :: imulti,get_mass

  if(get_mass)then
     write(csnap,'(I3.3)')isnap 
     write(logunit,*)'  -reading snapshot data (with masses): isnap = '//trim(adjustl(csnap))
  else
     write(csnap,'(I3.3)')isnap 
     write(logunit,*)'  -reading snapshot data (without masses): isnap = '//trim(adjustl(csnap))
  endif
  call flush(logunit)
  
  ! Look for snapshot as single file; if not there look for the .0 etc extension 
  ! Turn off hdf5 error messages
  call h5eset_auto_f(0,error)  
  write(suffix,'(I3.3)')isnap 
  infile_no_ext = trim(adjustl(basedir))//'/snapshot_'//trim(adjustl(suffix))&
       //'/'//trim(adjustl(basefile))//'_'//trim(adjustl(suffix))
  
  ! Test for presence of file - if not there, add on .0 to the filename
  snapfile = trim(adjustl(infile_no_ext))//'.hdf5'

  call h5fopen_f(snapfile,H5F_ACC_RDONLY_F,file_id,error)
  if(error.lt.0) then
     ! Can't open, try .0.hdf5 instead and try again
     snapfile = trim(adjustl(infile_no_ext))//'.0.hdf5'
     call h5fopen_f(snapfile,H5F_ACC_RDONLY_F,file_id,error)
     if(error < 0) then
        write(logunit,*)trim(snapfile)
        call abort('   -Cannot find snapshot in single- or multi-file format')
     else
        write(logunit,*)'   -opening snapshot in multi-file format'
        call flush(logunit)
        imulti = .TRUE.
        call h5fclose_f(file_id,error)
     endif
  else
     write(logunit,*)'   -opening snapshot in single-file format'
     call flush(logunit)
     imulti = .FALSE.
     call h5fclose_f(file_id,error)
  endif
  
  ! Turn on hdf5 error messages again
  call h5eset_auto_f(1,error)

  ! Read the header to get NumPart_ThisFile
  call read_header()  
  
  ! allocate particle data
  write(logunit,*)'   -snapshot contains:'
  do i=0,5
     write(logunit,*)'    -type = ',i,' NumPart = ',Header%NumPart_Total(i)
  enddo
  call flush(logunit)	   

  ! allocate memory
  call allocate_snapshot()
  
  ioffset(:) = 0
  do ifile = 1,Header%NumFilesPerSnapshot

     write(cfile,'(I3.3)')ifile
     if(ifile /= Header%NumFilesPerSnapshot)then
        write (logunit,'(A)',advance='no')char(13)//'    -read_snapshot: file = '//trim(adjustl(cfile))
        call flush(logunit)
     else
        write (logunit,'(A)',advance='yes')char(13)//'    -read_snapshot: file = '//trim(adjustl(cfile))
        call flush(logunit)
     endif

     if(imulti)then
        write(suffix,*)ifile-1
        snapfile = trim(adjustl(infile_no_ext))//'.'//trim(adjustl(suffix))//'.hdf5'
     else
        snapfile = trim(adjustl(infile_no_ext))//'.hdf5'
     endif
     
     call hdf5_open_file(file_handle,snapfile, readonly=.true.)
     call hdf5_read_attribute(file_handle,"/Header/NumPart_ThisFile",Header%NumPart_ThisFile)

     if(Header%NumPart_ThisFile(0) > 0)then

        varname = 'PartType0/ParticleIDs'
        call hdf5_read_data(file_handle,varname,Type0(ioffset(0):ioffset(0)+Header%NumPart_ThisFile(0)-1)%ParticleIDs)
        ! Two-dimensional arrays need special treatment
        allocate(real2d(3,0:Header%NumPart_ThisFile(0)-1))
        varname = 'PartType0/Coordinates'
        call hdf5_read_data(file_handle,varname,real2d)
        do i=0,Header%NumPart_ThisFile(0)-1
           Type0(ioffset(0)+i)%Coordinates(:) = real2d(:,i)
        enddo
        deallocate(real2d)

        if(get_mass)then
           if(Header%MassTable(0) > 0.0d0)then 
              Type0(ioffset(0):ioffset(0)+Header%NumPart_ThisFile(0)-1)%Mass = Header%MassTable(0)
           else
              varname = 'PartType0/Mass'
              call hdf5_read_data(file_handle,varname,Type0(ioffset(0):ioffset(0)+Header%NumPart_ThisFile(0)-1)%Mass)
           endif
        endif

     endif

     if(Header%NumPart_ThisFile(1) > 0)then

        varname = 'PartType1/ParticleIDs'
        call hdf5_read_data(file_handle,varname,Type1(ioffset(1):ioffset(1)+Header%NumPart_ThisFile(1)-1)%ParticleIDs)

        ! Two-dimensional arrays need special treatment
        allocate(real2d(3,0:Header%NumPart_ThisFile(1)-1))
        varname = 'PartType1/Coordinates'
        call hdf5_read_data(file_handle,varname,real2d)
        do i=0,Header%NumPart_ThisFile(1)-1
           Type1(ioffset(1)+i)%Coordinates(:) = real2d(:,i)
        enddo
        deallocate(real2d)

        if(get_mass)then
           if(Header%MassTable(1) > 0.0d0)then 
              Type1(ioffset(1):ioffset(1)+Header%NumPart_ThisFile(1)-1)%Mass = Header%MassTable(1)
           else
              varname = 'PartType1/Mass'
              call hdf5_read_data(file_handle,varname,Type1(ioffset(1):ioffset(1)+Header%NumPart_ThisFile(1)-1)%Mass)
           endif
        endif

     endif

     if(Header%NumPart_ThisFile(2) > 0)then

        varname = 'PartType2/ParticleIDs'
        call hdf5_read_data(file_handle,varname,Type2(ioffset(2):ioffset(2)+Header%NumPart_ThisFile(2)-1)%ParticleIDs)

        ! Two-dimensional arrays need special treatment
        allocate(real2d(3,0:Header%NumPart_ThisFile(2)-1))
        varname = 'PartType2/Coordinates'
        call hdf5_read_data(file_handle,varname,real2d)
        do i=0,Header%NumPart_ThisFile(2)-1
           Type2(ioffset(2)+i)%Coordinates(:) = real2d(:,i)
        enddo
        deallocate(real2d)

        if(get_mass)then
           if(Header%MassTable(2) > 0.0d0)then 
              Type2(ioffset(2):ioffset(2)+Header%NumPart_ThisFile(2)-1)%Mass = Header%MassTable(2)
           else
              varname = 'PartType2/Mass'
              call hdf5_read_data(file_handle,varname,Type2(ioffset(2):ioffset(2)+Header%NumPart_ThisFile(2)-1)%Mass)
           endif
        endif

     endif

     if(Header%NumPart_ThisFile(3) > 0)then

        varname = 'PartType3/ParticleIDs'
        call hdf5_read_data(file_handle,varname,Type3(ioffset(3):ioffset(3)+Header%NumPart_ThisFile(3)-1)%ParticleIDs)

        ! Two-dimensional arrays need special treatment
        allocate(real2d(3,0:Header%NumPart_ThisFile(3)-1))
        varname = 'PartType3/Coordinates'
        call hdf5_read_data(file_handle,varname,real2d)
        do i=0,Header%NumPart_ThisFile(3)-1
           Type3(ioffset(3)+i)%Coordinates(:) = real2d(:,i)
        enddo
        deallocate(real2d)

        if(get_mass)then
           if(Header%MassTable(3) > 0.0d0)then 
              Type3(ioffset(3):ioffset(3)+Header%NumPart_ThisFile(3)-1)%Mass = Header%MassTable(3)
           else
              varname = 'PartType3/Mass'
              call hdf5_read_data(file_handle,varname,Type3(ioffset(3):ioffset(3)+Header%NumPart_ThisFile(3)-1)%Mass)
           endif
        endif

     endif

     if(Header%NumPart_ThisFile(4) > 0)then

        varname = 'PartType4/ParticleIDs'
        call hdf5_read_data(file_handle,varname,Type4(ioffset(4):ioffset(4)+Header%NumPart_ThisFile(4)-1)%ParticleIDs)

        ! Two-dimensional arrays need special treatment
        allocate(real2d(3,0:Header%NumPart_ThisFile(4)-1))
        varname = 'PartType4/Coordinates'
        call hdf5_read_data(file_handle,varname,real2d)
        do i=0,Header%NumPart_ThisFile(4)-1
           Type4(ioffset(4)+i)%Coordinates(:) = real2d(:,i)
        enddo
        deallocate(real2d)

        if(get_mass)then
           if(Header%MassTable(4) > 0.0d0)then 
              Type3(ioffset(4):ioffset(4)+Header%NumPart_ThisFile(4)-1)%Mass = Header%MassTable(4)
           else
              varname = 'PartType4/Mass'
              call hdf5_read_data(file_handle,varname,Type4(ioffset(4):ioffset(4)+Header%NumPart_ThisFile(4)-1)%Mass)
           endif
        endif

     endif

     if(Header%NumPart_ThisFile(5) > 0)then

        varname = 'PartType5/ParticleIDs'
        call hdf5_read_data(file_handle,varname,Type5(ioffset(5):ioffset(5)+Header%NumPart_ThisFile(5)-1)%ParticleIDs)

        ! Two-dimensional arrays need special treatment
        allocate(real2d(3,0:Header%NumPart_ThisFile(5)-1))
        varname = 'PartType5/Coordinates'
        call hdf5_read_data(file_handle,varname,real2d)
        do i=0,Header%NumPart_ThisFile(5)-1
           Type5(ioffset(5)+i)%Coordinates(:) = real2d(:,i)
        enddo
        deallocate(real2d)

        if(get_mass)then
           if(Header%MassTable(5) > 0.0d0)then 
              Type3(ioffset(5):ioffset(5)+Header%NumPart_ThisFile(5)-1)%Mass = Header%MassTable(5)
           else
              varname = 'PartType5/Mass'
              call hdf5_read_data(file_handle,varname,Type5(ioffset(5):ioffset(5)+Header%NumPart_ThisFile(5)-1)%Mass)
           endif
        endif     

     endif

     ioffset(0) = ioffset(0) + Header%NumPart_ThisFile(0)
     ioffset(1) = ioffset(1) + Header%NumPart_ThisFile(1)
     ioffset(2) = ioffset(2) + Header%NumPart_ThisFile(2)
     ioffset(3) = ioffset(3) + Header%NumPart_ThisFile(3)
     ioffset(4) = ioffset(4) + Header%NumPart_ThisFile(4)
     ioffset(5) = ioffset(5) + Header%NumPart_ThisFile(5)
     call hdf5_close_file(file_handle)

  enddo

  write(logunit,*)'   -done'
  call flush(logunit)

end subroutine read_snapshot

! **********************************************************

subroutine dump_selection_ParticleIDs(isnap,is_final)
  use log
  use myfiles
  use particles
  use hdf5_wrapper

  implicit none
  integer                             :: file_handle
  integer*4                           :: itype,ipart,isnap
  real*4, allocatable, dimension(:,:) :: real2d
  character(len=512)                  :: PartTypeName
  logical, intent(in)                 :: is_final
  integer, dimension(0:5)             :: Num

  write(csnap,'(I3.3)')isnap
  if(is_final)then
     OutFile = trim(adjustl(OutDir))//'/OriginalSelection_'//trim(adjustl(csnap))//'.hdf5'
  else
     OutFile = trim(adjustl(OutDir))//'/SelectedParticles_'//trim(adjustl(csnap))//'.hdf5'
  endif

  write(logunit,*)' -dumping selection ParticleID...'
  call flush(logunit)

  call hdf5_create_file(file_handle,trim(adjustl(OutFile)))
  call hdf5_create_group(file_handle,'Header')
  call hdf5_write_attribute(file_handle,'Header/NumSelected',NumSelected)
  if(.NOT.(is_final)) call hdf5_write_attribute(file_handle,'Header/NumGathered',NumGathered)
  call hdf5_write_attribute(file_handle,'Header/x0',x0)
  call hdf5_write_attribute(file_handle,'Header/y0',y0)
  call hdf5_write_attribute(file_handle,'Header/z0',z0)
  call hdf5_write_attribute(file_handle,'Header/r0',r0)
  call hdf5_write_attribute(file_handle,'Header/Snapshot',isnap)

  if(is_final)then 
     Num = NumSelected
  else
     Num = NumGathered
  endif

  do itype=0,5
     if(Num(itype) > 0)then 
        selectcase(itype) 
        case(0)
           call hdf5_create_group(file_handle,'PartType0')
           call hdf5_write_data(file_handle,'PartType0/ParticleIDs',Type0_Selection(0:Num(0)-1)%ParticleIDs)
           call hdf5_write_data(file_handle,'PartType0/Index'      ,Type0_Selection(0:Num(0)-1)%Index)
             ! Two-dimensional arrays need special treatment
           allocate(real2d(3,0:Num(0)-1))
           do ipart=0,Num(0)-1
              real2d(:,ipart) = Type0_Selection(ipart)%Coordinates(:)
           enddo
           call hdf5_write_data(file_handle,'PartType0/Coordinates',real2d)
           deallocate(real2d)
        case(1)
           call hdf5_create_group(file_handle,'PartType1')
           call hdf5_write_data(file_handle,'PartType1/ParticleIDs',Type1_Selection(0:Num(1)-1)%ParticleIDs)
           call hdf5_write_data(file_handle,'PartType1/Index'      ,Type1_Selection(0:Num(1)-1)%Index)
             ! Two-dimensional arrays need special treatment
           allocate(real2d(3,0:Num(1)-1))
           do ipart=0,Num(1)-1
              real2d(:,ipart) = Type1_Selection(ipart)%Coordinates(:)
           enddo
           call hdf5_write_data(file_handle,'PartType1/Coordinates',real2d)
           deallocate(real2d)
        case(2)
           call hdf5_create_group(file_handle,'PartType2')
           call hdf5_write_data(file_handle,'PartType2/ParticleIDs',Type2_Selection(0:Num(2)-1)%ParticleIDs)
           call hdf5_write_data(file_handle,'PartType2/Index'      ,Type2_Selection(0:Num(2)-1)%Index)
             ! Two-dimensional arrays need special treatment
           allocate(real2d(3,0:Num(2)-1))
           do ipart=0,Num(2)-1
              real2d(:,ipart) = Type2_Selection(ipart)%Coordinates(:)
           enddo
           call hdf5_write_data(file_handle,'PartType2/Coordinates',real2d)
           deallocate(real2d)
        case(3)
           call hdf5_create_group(file_handle,'PartType3')
           call hdf5_write_data(file_handle,'PartType3/ParticleIDs',Type3_Selection(0:Num(3)-1)%ParticleIDs)
           call hdf5_write_data(file_handle,'PartType3/Index'      ,Type3_Selection(0:Num(3)-1)%Index)
             ! Two-dimensional arrays need special treatment
           allocate(real2d(3,0:Num(3)-1))
           do ipart=0,Num(3)-1
              real2d(:,ipart) = Type3_Selection(ipart)%Coordinates(:)
           enddo
           call hdf5_write_data(file_handle,'PartType3/Coordinates',real2d)
           deallocate(real2d)
        case(4)
           call hdf5_create_group(file_handle,'PartType4')
           call hdf5_write_data(file_handle,'PartType4/ParticleIDs',Type4_Selection(0:Num(4)-1)%ParticleIDs)
           call hdf5_write_data(file_handle,'PartType4/Index'      ,Type4_Selection(0:Num(4)-1)%Index)
             ! Two-dimensional arrays need special treatment
           allocate(real2d(3,0:Num(4)-1))
           do ipart=0,Num(4)-1
              real2d(:,ipart) = Type4_Selection(ipart)%Coordinates(:)
           enddo
           call hdf5_write_data(file_handle,'PartType4/Coordinates',real2d)
           deallocate(real2d)
        case(5)
           call hdf5_create_group(file_handle,'PartType5')
           call hdf5_write_data(file_handle,'PartType5/ParticleIDs',Type5_Selection(0:Num(5)-1)%ParticleIDs)
           call hdf5_write_data(file_handle,'PartType5/Index'      ,Type5_Selection(0:Num(5)-1)%Index)
             ! Two-dimensional arrays need special treatment
           allocate(real2d(3,0:Num(5)-1))
           do ipart=0,Num(5)-1
              real2d(:,ipart) = Type5_Selection(ipart)%Coordinates(:)
           enddo
           call hdf5_write_data(file_handle,'PartType5/Coordinates',real2d)
           deallocate(real2d)
        endselect
     endif
  enddo
  call hdf5_close_file(file_handle)

  write(logunit,*)'  -done'
  call flush(logunit)

end subroutine dump_selection_ParticleIDs

! **********************************************************

subroutine read_selection_ParticleIDs(isnap)
  use log
  use myfiles
  use particles
  use hdf5_wrapper

  implicit none
  integer   :: file_handle
  integer*4 :: itype,ipart,isnap
  real*4, allocatable, dimension(:,:) :: real2d
  character(len=512) :: PartTypeName

  write(csnap,'(I3.3)')isnap
  InFile = trim(adjustl(OutDir))//'/OriginalSelection_'//trim(adjustl(csnap))//'.hdf5'
    
  write(logunit,*)' -reading selection ParticleID...'//trim(InFile)
  call flush(logunit)

  call hdf5_open_file(file_handle,trim(adjustl(InFile)))
  call hdf5_read_attribute(file_handle,'Header/NumSelected',NumSelected)

  write(logunit,*)'  -type0:',NumSelected(0)
  write(logunit,*)'  -type1:',NumSelected(1)
  write(logunit,*)'  -type2:',NumSelected(2)
  write(logunit,*)'  -type3:',NumSelected(3)
  write(logunit,*)'  -type4:',NumSelected(4)
  write(logunit,*)'  -type5:',NumSelected(5)

  ! Allocate selection memory - don't re-use the other subroutine, as 
  ! when reading the data in, we already know how much memory is needed
  call allocate_selection_preknown()

  do itype=0,5
     if(NumSelected(itype) > 0)then 
        selectcase(itype) 
        case(0)
           call hdf5_read_data(file_handle,'PartType0/ParticleIDs',Type0_Selection(:)%ParticleIDs)
           call hdf5_read_data(file_handle,'PartType0/Index'      ,Type0_Selection(:)%Index)
             ! Two-dimensional arrays need special treatment
           allocate(real2d(3,0:NumSelected(0)-1))
           call hdf5_read_data(file_handle,'PartType0/Coordinates',real2d)
           do ipart=0,NumSelected(0)-1
              real2d(:,ipart) = Type0_Selection(ipart)%Coordinates(:)
           enddo
           deallocate(real2d)
        case(1)
           call hdf5_read_data(file_handle,'PartType1/ParticleIDs',Type1_Selection(:)%ParticleIDs)
           call hdf5_read_data(file_handle,'PartType1/Index'      ,Type1_Selection(:)%Index)
             ! Two-dimensional arrays need special treatment
           allocate(real2d(3,0:NumSelected(1)-1))
           call hdf5_read_data(file_handle,'PartType1/Coordinates',real2d)
           do ipart=0,NumSelected(1)-1
              real2d(:,ipart) = Type1_Selection(ipart)%Coordinates(:)
           enddo
           deallocate(real2d)        
        case(2)
           call hdf5_read_data(file_handle,'PartType2/ParticleIDs',Type2_Selection(:)%ParticleIDs)
           call hdf5_read_data(file_handle,'PartType2/Index'      ,Type2_Selection(:)%Index)
             ! Two-dimensional arrays need special treatment
           allocate(real2d(3,0:NumSelected(2)-1))
           call hdf5_read_data(file_handle,'PartType2/Coordinates',real2d)
           do ipart=0,NumSelected(2)-1
              real2d(:,ipart) = Type2_Selection(ipart)%Coordinates(:)
           enddo
           deallocate(real2d)
        case(3)
           call hdf5_read_data(file_handle,'PartType3/ParticleIDs',Type3_Selection(:)%ParticleIDs)
           call hdf5_read_data(file_handle,'PartType3/Index'      ,Type3_Selection(:)%Index)
             ! Two-dimensional arrays need special treatment
           allocate(real2d(3,0:NumSelected(3)-1))
           call hdf5_read_data(file_handle,'PartType3/Coordinates',real2d)
           do ipart=0,NumSelected(3)-1
              real2d(:,ipart) = Type3_Selection(ipart)%Coordinates(:)
           enddo
           deallocate(real2d)
        case(4)
           call hdf5_read_data(file_handle,'PartType4/ParticleIDs',Type4_Selection(:)%ParticleIDs)
           call hdf5_read_data(file_handle,'PartType4/Index'      ,Type4_Selection(:)%Index)
             ! Two-dimensional arrays need special treatment
           allocate(real2d(3,0:NumSelected(4)-1))
           call hdf5_read_data(file_handle,'PartType4/Coordinates',real2d)
           do ipart=0,NumSelected(4)-1
              real2d(:,ipart) = Type4_Selection(ipart)%Coordinates(:)
           enddo
           deallocate(real2d)
        case(5)
           call hdf5_read_data(file_handle,'PartType5/ParticleIDs',Type5_Selection(:)%ParticleIDs)
           call hdf5_read_data(file_handle,'PartType5/Index'      ,Type5_Selection(:)%Index)
             ! Two-dimensional arrays need special treatment
           allocate(real2d(3,0:NumSelected(5)-1))
           call hdf5_read_data(file_handle,'PartType5/Coordinates',real2d)
           do ipart=0,NumSelected(5)-1
              real2d(:,ipart) = Type5_Selection(ipart)%Coordinates(:)
           enddo
           deallocate(real2d)
        endselect
     endif
  enddo
  call hdf5_close_file(file_handle)

  write(logunit,*)'  -done'
  call flush(logunit)

end subroutine read_selection_ParticleIDs

! **********************************************************

subroutine dump_traced_com_and_indices(isnap)
   use log
   use myfiles
   use particles
   use hdf5_wrapper
 
   implicit none
   integer   :: file_handle
   integer*4 :: itype,ipart,isnap
   real*4, allocatable, dimension(:,:) :: real2d
   character(len=512) :: PartTypeName
     
 
   write(csnap,'(I3.3)')isnap
   OutFile = trim(adjustl(OutDir))//'/Traced_'//trim(adjustl(csnap))//'.hdf5'

   write(logunit,*)'  -dumping traced indices...'
   call flush(logunit)

   call hdf5_create_file(file_handle,trim(adjustl(OutFile)))
   call hdf5_create_group(file_handle,'Header')
   call hdf5_write_attribute(file_handle,'Header/NumSelected' ,NumSelected)
   call hdf5_write_attribute(file_handle,'Header/NumTraced'   ,NumTraced)
   call hdf5_write_attribute(file_handle,'Header/CenterOfMass',COM)
   call hdf5_write_attribute(file_handle,'Header/x0',x0)
   call hdf5_write_attribute(file_handle,'Header/y0',y0)
   call hdf5_write_attribute(file_handle,'Header/z0',z0)
   call hdf5_write_attribute(file_handle,'Header/r0',r0)

   do itype=0,5
      if(NumTraced(itype) > 0)then 
         selectcase(itype) 
         case(0)
            call hdf5_create_group(file_handle,'PartType0')
            call hdf5_write_data(file_handle,'PartType0/Index',Type0_Traced_Index(0:NumTraced(0)-1))
         case(1)
            call hdf5_create_group(file_handle,'PartType1')
            call hdf5_write_data(file_handle,'PartType1/Index',Type1_Traced_Index(0:NumTraced(1)-1))
         case(2)
            call hdf5_create_group(file_handle,'PartType2')
            call hdf5_write_data(file_handle,'PartType2/Index',Type2_Traced_Index(0:NumTraced(2)-1))
         case(3)
            call hdf5_create_group(file_handle,'PartType3')
            call hdf5_write_data(file_handle,'PartType3/Index',Type3_Traced_Index(0:NumTraced(3)-1))
         case(4)
            call hdf5_create_group(file_handle,'PartType4')
            call hdf5_write_data(file_handle,'PartType4/Index',Type4_Traced_Index(0:NumTraced(4)-1))
         case(5)
            call hdf5_create_group(file_handle,'PartType5')
            call hdf5_write_data(file_handle,'PartType5/Index',Type5_Traced_Index(0:NumTraced(5)-1))
         endselect
      endif
   enddo
   call hdf5_close_file(file_handle)

   write(logunit,*)'  -done'
   call flush(logunit)

 end subroutine dump_traced_com_and_indices

! **********************************************************

subroutine read_traced_com(isnap)
   use log
   use myfiles
   use particles
   use hdf5_wrapper
 
   implicit none
   integer   :: file_handle
   integer*4 :: itype,ipart,isnap
   real*4, allocatable, dimension(:,:) :: real2d
   character(len=512) :: PartTypeName
     
 
   write(csnap,'(I3.3)')isnap
   InFile = trim(adjustl(OutDir))//'/Traced_'//trim(adjustl(csnap))//'.hdf5'

   write(logunit,*)'  -reading traced centre of mass...'
   call flush(logunit)

   call hdf5_open_file(file_handle,trim(adjustl(InFile)))
   call hdf5_read_attribute(file_handle,'Header/NumSelected' ,NumSelected)
   call hdf5_read_attribute(file_handle,'Header/NumTraced'   ,NumTraced)
   call hdf5_read_attribute(file_handle,'Header/CenterOfMass',COM)
   call hdf5_close_file(file_handle)

   write(logunit,*)'  -done'
   call flush(logunit)

 end subroutine read_traced_com


! **********************************************************

 subroutine read_header
   use hdf5_wrapper
   use myfiles, only : snapfile
   use params_and_constants, only : Header
   implicit none
   integer*4 :: file_handle  
 
   call hdf5_open_file(file_handle,trim(snapfile),readonly=.true.)
   call hdf5_read_attribute(file_handle,'Header/RunLabel',header%RunLabel)
   call hdf5_read_attribute(file_handle,'Header/NumPart_ThisFile',header%NumPart_ThisFile)
   call hdf5_read_attribute(file_handle,'Header/NumPart_Total',header%NumPart_Total)
   call hdf5_read_attribute(file_handle,'Header/NumPart_Total_HighWord',header%NumPart_Total_HighWord)
   call hdf5_read_attribute(file_handle,'Header/MassTable',header%MassTable)
   call hdf5_read_attribute(file_handle,'Header/ExpansionFactor',header%ExpansionFactor)
   call hdf5_read_attribute(file_handle,'Header/Time_GYR',header%Time_GYR)
   call hdf5_read_attribute(file_handle,'Header/Redshift',header%Redshift)
   call hdf5_read_attribute(file_handle,'Header/BoxSize',header%BoxSize)
   call hdf5_read_attribute(file_handle,'Header/NumFilesPerSnapshot',header%NumFilesPerSnapshot)
   call hdf5_read_attribute(file_handle,'Header/Omega0',header%Omega0)
   call hdf5_read_attribute(file_handle,'Header/OmegaLambda',header%OmegaLambda)
   call hdf5_read_attribute(file_handle,'Header/OmegaBaryon',header%OmegaBaryon)
   call hdf5_read_attribute(file_handle,'Header/HubbleParam',header%HubbleParam)
   call hdf5_read_attribute(file_handle,'Header/Flag_Sfr',header%Flag_Sfr)
   call hdf5_read_attribute(file_handle,'Header/Flag_Cooling',header%Flag_Cooling)
   call hdf5_read_attribute(file_handle,'Header/Flag_StellarAge',header%Flag_StellarAge)
   call hdf5_read_attribute(file_handle,'Header/Flag_Metals',header%Flag_Metals)
   call hdf5_read_attribute(file_handle,'Header/Flag_Feedback',header%Flag_Feedback)
   call hdf5_close_file(file_handle)
 
 end subroutine read_header
 
! ************************************************

subroutine read_constants
  use hdf5_wrapper
  use myfiles, only : snapfile
  use params_and_constants, only : Constants
  implicit none
  integer*4 :: file_handle  


  call hdf5_open_file(file_handle,trim(snapfile),readonly=.true.)
  call hdf5_read_attribute(file_handle,'Constants/PI',Constants%PI)
  call hdf5_read_attribute(file_handle,'Constants/GAMMA',Constants%GAMMA)
  call hdf5_read_attribute(file_handle,'Constants/GRAVITY',Constants%GRAVITY)
  call hdf5_read_attribute(file_handle,'Constants/SOLAR_MASS',Constants%SOLAR_MASS)
  call hdf5_read_attribute(file_handle,'Constants/SOLAR_LUM',Constants%SOLAR_LUM)
  call hdf5_read_attribute(file_handle,'Constants/RAD_CONST',Constants%RAD_CONST)
  call hdf5_read_attribute(file_handle,'Constants/AVOGADRO',  Constants%AVOGADRO)
  call hdf5_read_attribute(file_handle,'Constants/BOLTZMANN',  Constants%BOLTZMANN)
  call hdf5_read_attribute(file_handle,'Constants/GAS_CONST',  Constants%GAS_CONST)
  call hdf5_read_attribute(file_handle,'Constants/C', Constants%C)
  call hdf5_read_attribute(file_handle,'Constants/PLANCK',Constants%PLANCK)
  call hdf5_read_attribute(file_handle,'Constants/CM_PER_MPC',Constants%CM_PER_MPC)
  call hdf5_read_attribute(file_handle,'Constants/PROTONMASS',Constants%PROTONMASS)
  call hdf5_read_attribute(file_handle,'Constants/ELECTRONMASS',Constants%ELECTRONMASS)
  call hdf5_read_attribute(file_handle,'Constants/ELECTRONCHARGE',Constants%ELECTRONCHARGE)
  call hdf5_read_attribute(file_handle,'Constants/HUBBLE',Constants%HUBBLE)
  call hdf5_read_attribute(file_handle,'Constants/T_CMB0',Constants%T_CMB0)
  call hdf5_read_attribute(file_handle,'Constants/SEC_PER_MEGAYEAR',Constants%SEC_PER_MEGAYEAR)
  call hdf5_read_attribute(file_handle,'Constants/SEC_PER_YEAR',Constants%SEC_PER_YEAR)
  call hdf5_read_attribute(file_handle,'Constants/STEFAN',Constants%STEFAN)
  call hdf5_read_attribute(file_handle,'Constants/THOMPSON',Constants%THOMPSON)
  call hdf5_read_attribute(file_handle,'Constants/EV_TO_ERG',Constants%EV_TO_ERG)
  call hdf5_read_attribute(file_handle,'Constants/Z_Solar',Constants%Z_Solar)
  call hdf5_close_file(file_handle)

end subroutine read_constants

! ************************************************ 
 
 subroutine read_units
   use hdf5_wrapper
   use myfiles, only : snapfile
   use params_and_constants, only : Units
   implicit none
   integer*4 :: file_handle  
 
   call hdf5_open_file(file_handle,trim(snapfile),readonly=.true.)
   call hdf5_read_attribute(file_handle,'Units/UnitLength_in_cm',Units%UnitLength_in_cm)
   call hdf5_read_attribute(file_handle,'Units/UnitMass_in_g',Units%UnitMass_in_g)
   call hdf5_read_attribute(file_handle,'Units/UnitVelocity_in_cm_per_s',Units%UnitVelocity_in_cm_per_s)
   call hdf5_read_attribute(file_handle,'Units/UnitDensity_in_cgs',Units%UnitDensity_in_cgs)
   call hdf5_read_attribute(file_handle,'Units/UnitEnergy_in_cgs',Units%UnitEnergy_in_cgs)
   call hdf5_read_attribute(file_handle,'Units/UnitPressure_in_cgs',Units%UnitPressure_in_cgs)
   call hdf5_read_attribute(file_handle,'Units/UnitTime_in_s',Units%UnitTime_in_s)
   call hdf5_close_file(file_handle)
 
 end subroutine read_units
 
 ! ************************************************

 subroutine read_parameters
  use hdf5_wrapper
  use myfiles, only : snapfile
  use params_and_constants, only : Parameters
  implicit none
  integer*4 :: file_handle  

  call hdf5_open_file(file_handle,trim(snapfile),readonly=.true.)
  ! Chemical Elements
  call hdf5_read_attribute(file_handle,'Parameters/ChemicalElements/BG_NELEMENTS',&
       Parameters%BG_NELEMENTS)
  call hdf5_read_attribute(file_handle,'Parameters/ChemicalElements/ElementNames',&
       Parameters%ElementNames)
  call hdf5_read_attribute(file_handle,'Parameters/ChemicalElements/InitAbundance_Hydrogen',&
       Parameters%InitAbundance_Hydrogen)
  call hdf5_read_attribute(file_handle,'Parameters/ChemicalElements/InitAbundance_Helium',&
       Parameters%InitAbundance_Helium)
  call hdf5_read_attribute(file_handle,'Parameters/ChemicalElements/InitAbundance_Carbon',&
       Parameters%InitAbundance_Carbon)
  call hdf5_read_attribute(file_handle,'Parameters/ChemicalElements/InitAbundance_Nitrogen',&
       Parameters%InitAbundance_Nitrogen)
  call hdf5_read_attribute(file_handle,'Parameters/ChemicalElements/InitAbundance_Oxygen',&
       Parameters%InitAbundance_Oxygen)
  call hdf5_read_attribute(file_handle,'Parameters/ChemicalElements/InitAbundance_Neon',&
       Parameters%InitAbundance_Neon)
  call hdf5_read_attribute(file_handle,'Parameters/ChemicalElements/InitAbundance_Magnesium',&
       Parameters%InitAbundance_Magnesium)
  call hdf5_read_attribute(file_handle,'Parameters/ChemicalElements/InitAbundance_Silicon',&
       Parameters%InitAbundance_Silicon)
  call hdf5_read_attribute(file_handle,'Parameters/ChemicalElements/InitAbundance_Iron',&
       Parameters%InitAbundance_Iron)
  call hdf5_read_attribute(file_handle,'Parameters/ChemicalElements/CalciumOverSilicon',&
       Parameters%CalciumOverSilicon)
  call hdf5_read_attribute(file_handle,'Parameters/ChemicalElements/SulphurOverSilicon',&
       Parameters%SulphurOverSilicon)

!  ! Numerical Parameters
!  call hdf5_read_attribute(file_handle,'Parameters/NumericalParameters/ComovingIntegrationOn',&
!       Parameters%ComovingIntegrationOn)
!  call hdf5_read_attribute(file_handle,'Parameters/NumericalParameters/TypeOfTimestepCriterion',&
!       Parameters%TypeOfTimeStepCriterion)
!  call hdf5_read_attribute(file_handle,'Parameters/NumericalParameters/ErrTolIntAccuracy',&
!       Parameters%ErrTolIntAccuracy)
!  call hdf5_read_attribute(file_handle,'Parameters/NumericalParameters/MaxSizeTimestep',&
!       Parameters%MaxSizeTimestep)
!  call hdf5_read_attribute(file_handle,'Parameters/NumericalParameters/MinSizeTimestep',&
!       Parameters%MinSizeTimestep)
!  call hdf5_read_attribute(file_handle,'Parameters/NumericalParameters/ErrTolTheta',&
!       Parameters%ErrTolTheta)
!  call hdf5_read_attribute(file_handle,'Parameters/NumericalParameters/TypeOfOpeningCriterion',&
!       Parameters%TypeOfOpeningCriterion)
!  call hdf5_read_attribute(file_handle,'Parameters/NumericalParameters/ErrTolForceAcc',&
!       Parameters%ErrTolForceAcc)
!  call hdf5_read_attribute(file_handle,'Parameters/NumericalParameters/TreeDomainUpdateFrequency',&
!       Parameters%TreeDomainUpdateFrequency)
!  call hdf5_read_attribute(file_handle,'Parameters/NumericalParameters/DesNumNgb',&
!       Parameters%DesNumNgb)
!  call hdf5_read_attribute(file_handle,'Parameters/NumericalParameters/MaxNumNgbDeviation',&
!       Parameters%MaxNumNgbDeviation)
!  call hdf5_read_attribute(file_handle,'Parameters/NumericalParameters/ArtBulkViscConst',&
!       Parameters%ArtBulkViscConst)
!  call hdf5_read_attribute(file_handle,'Parameters/NumericalParameters/InitGasU_ERG',&
!       Parameters%InitGasU_ERG)
!  call hdf5_read_attribute(file_handle,'Parameters/NumericalParameters/MinGasU_ERG',&
!       Parameters%MinGasU_ERG)
!  call hdf5_read_attribute(file_handle,'Parameters/NumericalParameters/CourantFac',&
!       Parameters%CourantFac)
!  call hdf5_read_attribute(file_handle,'Parameters/NumericalParameters/PartAllocFactor',&
!       Parameters%PartAllocFactor)
!  call hdf5_read_attribute(file_handle,'Parameters/NumericalParameters/TreeAllocFactor',&
!       Parameters%TreeAllocFactor)
!  call hdf5_read_attribute(file_handle,'Parameters/NumericalParameters/BufferSize',&
!       Parameters%BufferSize)
!  call hdf5_read_attribute(file_handle,'Parameters/NumericalParameters/MinGasHsmlFractional',&
!       Parameters%MinGasHsmlFractional)
!  
!  ! Stellar Evolution Parameters
!  call hdf5_read_attribute(file_handle,'Parameters/StellarEvolutionParameters/SF_EOSGammaEffective',&
!       Parameters%SF_EOSGammaEffective)
!  call hdf5_read_attribute(file_handle,'Parameters/StellarEvolutionParameters/SF_EOSEnergyAtThreshold_ERG',&
!       Parameters%SF_EOSEnergyAtThreshold_ERG)
!  call hdf5_read_attribute(file_handle,'Parameters/StellarEvolutionParameters/SF_THRESH_MinPhysDens_HpCM3',&
!       Parameters%SF_THRESH_MinPhysDens_HpCM3)
!  call hdf5_read_attribute(file_handle,'Parameters/StellarEvolutionParameters/SF_THRESH_MinOverDens',&
!       Parameters%SF_THRESH_MinOverDens)
!
!  ! Seems to have gone in new code
!  !call hdf5_read_attribute(file_handle,'Parameters/StellarEvolutionParameters/SF_THRESH_MetDepExponent',&
!  !     Parameters%SF_THRESH_MetDepExponent)
!  !call hdf5_read_attribute(file_handle,'Parameters/StellarEvolutionParameters/SF_THRESH_MetDepFiducialZ',&
!  !     Parameters%SF_THRESH_MetDepFiducialZ)
!  !call hdf5_read_attribute(file_handle,'Parameters/StellarEvolutionParameters/SF_THRESH_MetDepMaxPhysDens_HpCM3',&
!  !     Parameters%SF_THRESH_MetDepMaxPhysDens_HpCM3)
!  call hdf5_read_attribute(file_handle,'Parameters/StellarEvolutionParameters/SF_THRESH_MaxTemp_K',&
!       Parameters%SF_THRESH_MaxTemp_K)
!
!  call hdf5_read_attribute(file_handle,'Parameters/StellarEvolutionParameters/SF_SchmidtLawCoeff_MSUNpYRpKPC2',&
!       Parameters%SF_SchmidtLawCoeff_MSUNpYRpKPC2)
!  call hdf5_read_attribute(file_handle,'Parameters/StellarEvolutionParameters/SF_SchmidtLawExponent',&
!       Parameters%SF_SchmidtLawExponent)
!  call hdf5_read_attribute(file_handle,'Parameters/StellarEvolutionParameters/SF_SchmidtLawCoeff_GpSpCM2',&
!       Parameters%SF_SchmidtLawCoeff_GpSpCM2)
!  call hdf5_read_attribute(file_handle,'Parameters/StellarEvolutionParameters/SF_MetDepDensThresholdOn',&
!       Parameters%SF_MetDepDensThresholdOn)
!  call hdf5_read_attribute(file_handle,'Parameters/StellarEvolutionParameters/IMF_Model',&
!       Parameters%IMF_Model)
!  call hdf5_read_attribute(file_handle,'Parameters/StellarEvolutionParameters/IMF_LifetimeModel',&
!       Parameters%IMF_LifetimeModel)
!  call hdf5_read_attribute(file_handle,'Parameters/StellarEvolutionParameters/IMF_MinMass_MSUN',&
!       Parameters%IMF_MinMass_MSUN)
!  call hdf5_read_attribute(file_handle,'Parameters/StellarEvolutionParameters/IMF_MaxMass_MSUN',&
!       Parameters%IMF_MaxMass_MSUN)
!  call hdf5_read_attribute(file_handle,'Parameters/StellarEvolutionParameters/SNIa_Model',&
!       Parameters%SNIa_Model)
!  call hdf5_read_attribute(file_handle,'Parameters/StellarEvolutionParameters/SNIa_Efficiency_fracwd',&
!       Parameters%SNIa_Efficiency_fracwd)
!  call hdf5_read_attribute(file_handle,'Parameters/StellarEvolutionParameters/SNIa_Energy_ERG',&
!       Parameters%SNIa_Energy_ERG)
!  call hdf5_read_attribute(file_handle,'Parameters/StellarEvolutionParameters/SNIa_MassTransferOn',&
!       Parameters%SNIa_MassTransferOn)
!  call hdf5_read_attribute(file_handle,'Parameters/StellarEvolutionParameters/SNIa_EnergyTransferOn',&
!       Parameters%SNIa_EnergyTransferOn)
!  call hdf5_read_attribute(file_handle,'Parameters/StellarEvolutionParameters/SNII_MassTransferOn',&
!       Parameters%SNII_MassTransferOn)
!  call hdf5_read_attribute(file_handle,'Parameters/StellarEvolutionParameters/SNII_MinMass_MSUN',&
!       Parameters%SNII_MinMass_MSUN)
!  call hdf5_read_attribute(file_handle,'Parameters/StellarEvolutionParameters/SNII_MaxMass_MSUN',&
!       Parameters%SNII_MaxMass_MSUN)
!  call hdf5_read_attribute(file_handle,'Parameters/StellarEvolutionParameters/SNII_Factor_Hydrogen',&
!       Parameters%SNII_Factor_Hydrogen)
!  call hdf5_read_attribute(file_handle,'Parameters/StellarEvolutionParameters/SNII_Factor_Helium',&
!       Parameters%SNII_Factor_Helium)
!  call hdf5_read_attribute(file_handle,'Parameters/StellarEvolutionParameters/SNII_Factor_Carbon',&
!       Parameters%SNII_Factor_Carbon)
!  call hdf5_read_attribute(file_handle,'Parameters/StellarEvolutionParameters/SNII_Factor_Nitrogen',&
!       Parameters%SNII_Factor_Nitrogen)
!  call hdf5_read_attribute(file_handle,'Parameters/StellarEvolutionParameters/SNII_Factor_Oxygen',&
!       Parameters%SNII_Factor_Oxygen)
!  call hdf5_read_attribute(file_handle,'Parameters/StellarEvolutionParameters/SNII_Factor_Neon',&
!       Parameters%SNII_Factor_Neon)
!  call hdf5_read_attribute(file_handle,'Parameters/StellarEvolutionParameters/SNII_Factor_Magnesium',&
!       Parameters%SNII_Factor_Magnesium)
!  call hdf5_read_attribute(file_handle,'Parameters/StellarEvolutionParameters/SNII_Factor_Silicon',&
!       Parameters%SNII_Factor_Silicon)
!  call hdf5_read_attribute(file_handle,'Parameters/StellarEvolutionParameters/SNII_Factor_Iron',&
!       Parameters%SNII_Factor_Iron)
!  call hdf5_read_attribute(file_handle,'Parameters/StellarEvolutionParameters/AGB_MassTransferOn',&
!       Parameters%AGB_MassTransferOn)
!  call hdf5_read_attribute(file_handle,'Parameters/StellarEvolutionParameters/POPIII_MassTransferOn',&
!       Parameters%POPIII_MassTransferOn)
!  call hdf5_read_attribute(file_handle,'Parameters/StellarEvolutionParameters/POPIII_EnergyTransferOn',&
!       Parameters%POPIII_EnergyTransferOn)
!  call hdf5_read_attribute(file_handle,'Parameters/StellarEvolutionParameters/POPIII_Energy_ERG',&
!       Parameters%POPIII_Energy_ERG)
!  call hdf5_read_attribute(file_handle,'Parameters/StellarEvolutionParameters/POPIII_NumPerMsun',&
!       Parameters%POPIII_NumPerMsun)
!
!  ! Wind Parameters
!  call hdf5_read_attribute(file_handle,'Parameters/WindParameters/SNII_WindOn',&
!       Parameters%SNII_WindOn)
!  call hdf5_read_attribute(file_handle,'Parameters/WindParameters/SNII_WindIsotropicOn',&
!       Parameters%SNII_WindIsotropicOn)
!  call hdf5_read_attribute(file_handle,'Parameters/WindParameters/SNII_WindSpeed_KMpS',&
!       Parameters%SNII_WindSpeed_KMpS)
!  call hdf5_read_attribute(file_handle,'Parameters/WindParameters/SNII_WindMassLoading',&
!       Parameters%SNII_WindMassLoading)
!  call hdf5_read_attribute(file_handle,'Parameters/WindParameters/SNII_WindDelay_YR',&
!       Parameters%SNII_WindDelay_YR)

  call hdf5_close_file(file_handle)

end subroutine read_parameters


! ************************************************
!!! 
!!! subroutine write_parameters(outfile)
!!!   use hdf5_wrapper
!!!   use params_and_constants, only : Parameters
!!!   implicit none
!!!   integer*4 :: file_handle  
!!!   character(len=*) :: outfile
!!! 
!!!   call hdf5_open_file(file_handle,trim(outfile))
!!! 
!!!   ! Create groups
!!!   call hdf5_create_group(file_handle,'Parameters')
!!!   call hdf5_create_group(file_handle,'Parameters/ChemicalElements')
!!!   call hdf5_create_group(file_handle,'Parameters/NumericalParameters')
!!!   call hdf5_create_group(file_handle,'Parameters/StellarEvolutionParameters')
!!!   call hdf5_create_group(file_handle,'Parameters/WindParameters')
!!! 
!!!   ! Chemical Elements
!!!   call hdf5_write_attribute(file_handle,'Parameters/ChemicalElements/BG_NELEMENTS',&
!!!        Parameters%BG_NELEMENTS)
!!!   call hdf5_write_attribute(file_handle,'Parameters/ChemicalElements/ElementNames',&
!!!        Parameters%ElementNames)
!!!   call hdf5_write_attribute(file_handle,'Parameters/ChemicalElements/InitAbundance_Hydrogen',&
!!!        Parameters%InitAbundance_Hydrogen)
!!!   call hdf5_write_attribute(file_handle,'Parameters/ChemicalElements/InitAbundance_Helium',&
!!!        Parameters%InitAbundance_Helium)
!!!   call hdf5_write_attribute(file_handle,'Parameters/ChemicalElements/InitAbundance_Carbon',&
!!!        Parameters%InitAbundance_Carbon)
!!!   call hdf5_write_attribute(file_handle,'Parameters/ChemicalElements/InitAbundance_Nitrogen',&
!!!        Parameters%InitAbundance_Nitrogen)
!!!   call hdf5_write_attribute(file_handle,'Parameters/ChemicalElements/InitAbundance_Oxygen',&
!!!        Parameters%InitAbundance_Oxygen)
!!!   call hdf5_write_attribute(file_handle,'Parameters/ChemicalElements/InitAbundance_Neon',&
!!!        Parameters%InitAbundance_Neon)
!!!   call hdf5_write_attribute(file_handle,'Parameters/ChemicalElements/InitAbundance_Magnesium',&
!!!        Parameters%InitAbundance_Magnesium)
!!!   call hdf5_write_attribute(file_handle,'Parameters/ChemicalElements/InitAbundance_Silicon',&
!!!        Parameters%InitAbundance_Silicon)
!!!   call hdf5_write_attribute(file_handle,'Parameters/ChemicalElements/InitAbundance_Iron',&
!!!        Parameters%InitAbundance_Iron)
!!!   call hdf5_write_attribute(file_handle,'Parameters/ChemicalElements/CalciumOverSilicon',&
!!!        Parameters%CalciumOverSilicon)
!!!   call hdf5_write_attribute(file_handle,'Parameters/ChemicalElements/SulphurOverSilicon',&
!!!        Parameters%SulphurOverSilicon)
!!! 
!!!   ! Numerical Parameters
!!!   call hdf5_write_attribute(file_handle,'Parameters/NumericalParameters/ComovingIntegrationOn',&
!!!        Parameters%ComovingIntegrationOn)
!!!   call hdf5_write_attribute(file_handle,'Parameters/NumericalParameters/TypeOfTimestepCriterion',&
!!!        Parameters%TypeOfTimeStepCriterion)
!!!   call hdf5_write_attribute(file_handle,'Parameters/NumericalParameters/ErrTolIntAccuracy',&
!!!        Parameters%ErrTolIntAccuracy)
!!!   call hdf5_write_attribute(file_handle,'Parameters/NumericalParameters/MaxSizeTimestep',&
!!!        Parameters%MaxSizeTimestep)
!!!   call hdf5_write_attribute(file_handle,'Parameters/NumericalParameters/MinSizeTimestep',&
!!!        Parameters%MinSizeTimestep)
!!!   call hdf5_write_attribute(file_handle,'Parameters/NumericalParameters/ErrTolTheta',&
!!!        Parameters%ErrTolTheta)
!!!   call hdf5_write_attribute(file_handle,'Parameters/NumericalParameters/TypeOfOpeningCriterion',&
!!!        Parameters%TypeOfOpeningCriterion)
!!!   call hdf5_write_attribute(file_handle,'Parameters/NumericalParameters/ErrTolForceAcc',&
!!!        Parameters%ErrTolForceAcc)
!!!   call hdf5_write_attribute(file_handle,'Parameters/NumericalParameters/TreeDomainUpdateFrequency',&
!!!        Parameters%TreeDomainUpdateFrequency)
!!!   call hdf5_write_attribute(file_handle,'Parameters/NumericalParameters/DesNumNgb',&
!!!        Parameters%DesNumNgb)
!!!   call hdf5_write_attribute(file_handle,'Parameters/NumericalParameters/MaxNumNgbDeviation',&
!!!        Parameters%MaxNumNgbDeviation)
!!!   call hdf5_write_attribute(file_handle,'Parameters/NumericalParameters/ArtBulkViscConst',&
!!!        Parameters%ArtBulkViscConst)
!!!   call hdf5_write_attribute(file_handle,'Parameters/NumericalParameters/InitGasU_ERG',&
!!!        Parameters%InitGasU_ERG)
!!!   call hdf5_write_attribute(file_handle,'Parameters/NumericalParameters/MinGasU_ERG',&
!!!        Parameters%MinGasU_ERG)
!!!   call hdf5_write_attribute(file_handle,'Parameters/NumericalParameters/CourantFac',&
!!!        Parameters%CourantFac)
!!!   call hdf5_write_attribute(file_handle,'Parameters/NumericalParameters/PartAllocFactor',&
!!!        Parameters%PartAllocFactor)
!!!   call hdf5_write_attribute(file_handle,'Parameters/NumericalParameters/TreeAllocFactor',&
!!!        Parameters%TreeAllocFactor)
!!!   call hdf5_write_attribute(file_handle,'Parameters/NumericalParameters/BufferSize',&
!!!        Parameters%BufferSize)
!!!   call hdf5_write_attribute(file_handle,'Parameters/NumericalParameters/MinGasHsmlFractional',&
!!!        Parameters%MinGasHsmlFractional)
!!! 
!!!   ! Stellar Evolution Parameters
!!!   call hdf5_write_attribute(file_handle,'Parameters/StellarEvolutionParameters/SF_EOSGammaEffective',&
!!!        Parameters%SF_EOSGammaEffective)
!!!   call hdf5_write_attribute(file_handle,'Parameters/StellarEvolutionParameters/SF_EOSEnergyAtThreshold_ERG',&
!!!        Parameters%SF_EOSEnergyAtThreshold_ERG)
!!!   call hdf5_write_attribute(file_handle,'Parameters/StellarEvolutionParameters/SF_THRESH_MinPhysDens_HpCM3',&
!!!        Parameters%SF_THRESH_MinPhysDens_HpCM3)
!!!   call hdf5_write_attribute(file_handle,'Parameters/StellarEvolutionParameters/SF_THRESH_MinOverDens',&
!!!        Parameters%SF_THRESH_MinOverDens)
!!!   call hdf5_write_attribute(file_handle,'Parameters/StellarEvolutionParameters/SF_THRESH_MetDepExponent',&
!!!        Parameters%SF_THRESH_MetDepExponent)
!!!   call hdf5_write_attribute(file_handle,'Parameters/StellarEvolutionParameters/SF_THRESH_MetDepFiducialZ',&
!!!        Parameters%SF_THRESH_MetDepFiducialZ)
!!!   call hdf5_write_attribute(file_handle,'Parameters/StellarEvolutionParameters/SF_THRESH_MetDepMaxPhysDens_HpCM3',&
!!!        Parameters%SF_THRESH_MetDepMaxPhysDens_HpCM3)
!!!   call hdf5_write_attribute(file_handle,'Parameters/StellarEvolutionParameters/SF_THRESH_MaxTemp_K',&
!!!        Parameters%SF_THRESH_MaxTemp_K)
!!!   call hdf5_write_attribute(file_handle,'Parameters/StellarEvolutionParameters/SF_SchmidtLawCoeff_MSUNpYRpKPC2',&
!!!        Parameters%SF_SchmidtLawCoeff_MSUNpYRpKPC2)
!!!   call hdf5_write_attribute(file_handle,'Parameters/StellarEvolutionParameters/SF_SchmidtLawExponent',&
!!!        Parameters%SF_SchmidtLawExponent)
!!!   call hdf5_write_attribute(file_handle,'Parameters/StellarEvolutionParameters/SF_SchmidtLawCoeff_GpSpCM2',&
!!!        Parameters%SF_SchmidtLawCoeff_GpSpCM2)
!!!   call hdf5_write_attribute(file_handle,'Parameters/StellarEvolutionParameters/SF_MetDepDensThresholdOn',&
!!!        Parameters%SF_MetDepDensThresholdOn)
!!!   call hdf5_write_attribute(file_handle,'Parameters/StellarEvolutionParameters/IMF_Model',&
!!!        Parameters%IMF_Model)
!!!   call hdf5_write_attribute(file_handle,'Parameters/StellarEvolutionParameters/IMF_LifetimeModel',&
!!!        Parameters%IMF_LifetimeModel)
!!!   call hdf5_write_attribute(file_handle,'Parameters/StellarEvolutionParameters/IMF_MinMass_MSUN',&
!!!        Parameters%IMF_MinMass_MSUN)
!!!   call hdf5_write_attribute(file_handle,'Parameters/StellarEvolutionParameters/IMF_MaxMass_MSUN',&
!!!        Parameters%IMF_MaxMass_MSUN)
!!!   call hdf5_write_attribute(file_handle,'Parameters/StellarEvolutionParameters/SNIa_Model',&
!!!        Parameters%SNIa_Model)
!!!   call hdf5_write_attribute(file_handle,'Parameters/StellarEvolutionParameters/SNIa_Efficiency_fracwd',&
!!!        Parameters%SNIa_Efficiency_fracwd)
!!!   call hdf5_write_attribute(file_handle,'Parameters/StellarEvolutionParameters/SNIa_Energy_ERG',&
!!!        Parameters%SNIa_Energy_ERG)
!!!   call hdf5_write_attribute(file_handle,'Parameters/StellarEvolutionParameters/SNIa_MassTransferOn',&
!!!        Parameters%SNIa_MassTransferOn)
!!!   call hdf5_write_attribute(file_handle,'Parameters/StellarEvolutionParameters/SNIa_EnergyTransferOn',&
!!!        Parameters%SNIa_EnergyTransferOn)
!!!   call hdf5_write_attribute(file_handle,'Parameters/StellarEvolutionParameters/SNII_MassTransferOn',&
!!!        Parameters%SNII_MassTransferOn)
!!!   call hdf5_write_attribute(file_handle,'Parameters/StellarEvolutionParameters/SNII_MinMass_MSUN',&
!!!        Parameters%SNII_MinMass_MSUN)
!!!   call hdf5_write_attribute(file_handle,'Parameters/StellarEvolutionParameters/SNII_MaxMass_MSUN',&
!!!        Parameters%SNII_MaxMass_MSUN)
!!!   call hdf5_write_attribute(file_handle,'Parameters/StellarEvolutionParameters/SNII_Factor_Hydrogen',&
!!!        Parameters%SNII_Factor_Hydrogen)
!!!   call hdf5_write_attribute(file_handle,'Parameters/StellarEvolutionParameters/SNII_Factor_Helium',&
!!!        Parameters%SNII_Factor_Helium)
!!!   call hdf5_write_attribute(file_handle,'Parameters/StellarEvolutionParameters/SNII_Factor_Carbon',&
!!!        Parameters%SNII_Factor_Carbon)
!!!   call hdf5_write_attribute(file_handle,'Parameters/StellarEvolutionParameters/SNII_Factor_Nitrogen',&
!!!        Parameters%SNII_Factor_Nitrogen)
!!!   call hdf5_write_attribute(file_handle,'Parameters/StellarEvolutionParameters/SNII_Factor_Oxygen',&
!!!        Parameters%SNII_Factor_Oxygen)
!!!   call hdf5_write_attribute(file_handle,'Parameters/StellarEvolutionParameters/SNII_Factor_Neon',&
!!!        Parameters%SNII_Factor_Neon)
!!!   call hdf5_write_attribute(file_handle,'Parameters/StellarEvolutionParameters/SNII_Factor_Magnesium',&
!!!        Parameters%SNII_Factor_Magnesium)
!!!   call hdf5_write_attribute(file_handle,'Parameters/StellarEvolutionParameters/SNII_Factor_Silicon',&
!!!        Parameters%SNII_Factor_Silicon)
!!!   call hdf5_write_attribute(file_handle,'Parameters/StellarEvolutionParameters/SNII_Factor_Iron',&
!!!        Parameters%SNII_Factor_Iron)
!!!   call hdf5_write_attribute(file_handle,'Parameters/StellarEvolutionParameters/AGB_MassTransferOn',&
!!!        Parameters%AGB_MassTransferOn)
!!!   call hdf5_write_attribute(file_handle,'Parameters/StellarEvolutionParameters/POPIII_MassTransferOn',&
!!!        Parameters%POPIII_MassTransferOn)
!!!   call hdf5_write_attribute(file_handle,'Parameters/StellarEvolutionParameters/POPIII_EnergyTransferOn',&
!!!        Parameters%POPIII_EnergyTransferOn)
!!!   call hdf5_write_attribute(file_handle,'Parameters/StellarEvolutionParameters/POPIII_Energy_ERG',&
!!!        Parameters%POPIII_Energy_ERG)
!!!   call hdf5_write_attribute(file_handle,'Parameters/StellarEvolutionParameters/POPIII_NumPerMsun',&
!!!        Parameters%POPIII_NumPerMsun)
!!! 
!!!   ! Wind Parameters
!!!   call hdf5_write_attribute(file_handle,'Parameters/WindParameters/SNII_WindOn',&
!!!        Parameters%SNII_WindOn)
!!!   call hdf5_write_attribute(file_handle,'Parameters/WindParameters/SNII_WindIsotropicOn',&
!!!        Parameters%SNII_WindIsotropicOn)
!!!   call hdf5_write_attribute(file_handle,'Parameters/WindParameters/SNII_WindSpeed_KMpS',&
!!!        Parameters%SNII_WindSpeed_KMpS)
!!!   call hdf5_write_attribute(file_handle,'Parameters/WindParameters/SNII_WindMassLoading',&
!!!        Parameters%SNII_WindMassLoading)
!!!   call hdf5_write_attribute(file_handle,'Parameters/WindParameters/SNII_WindDelay_YR',&
!!!        Parameters%SNII_WindDelay_YR)
!!! 
!!!   call hdf5_close_file(file_handle)
!!! 
!!! end subroutine write_parameters
!!! 
!!! 
! ************************************************

subroutine read_cgs_conversions
  use hdf5_wrapper
  use myfiles, only : snapfile
  use params_and_constants, only : Conv
  implicit none
  integer*4 :: file_handle  

  call hdf5_open_file(file_handle,trim(snapfile),readonly=.true.)
  call hdf5_read_attribute(file_handle,'PartType0/Coordinates/CGSConversionFactor',Conv%Coordinates%CGSConversionFactor)
  call hdf5_read_attribute(file_handle,'PartType0/Coordinates/aexp-scale-exponent',Conv%Coordinates%aexp_scale_exponent)
  call hdf5_read_attribute(file_handle,'PartType0/Coordinates/h-scale-exponent'   ,Conv%Coordinates%h_scale_exponent)
  call hdf5_read_attribute(file_handle,'PartType0/Velocity/CGSConversionFactor'   ,Conv%Velocity%CGSConversionFactor)
  call hdf5_read_attribute(file_handle,'PartType0/Velocity/aexp-scale-exponent'   ,Conv%Velocity%aexp_scale_exponent)
  call hdf5_read_attribute(file_handle,'PartType0/Velocity/h-scale-exponent'      ,Conv%Velocity%h_scale_exponent)
  call hdf5_read_attribute(file_handle,'PartType0/Temperature/CGSConversionFactor',Conv%Temperature%CGSConversionFactor)
  call hdf5_read_attribute(file_handle,'PartType0/Temperature/aexp-scale-exponent',Conv%Temperature%aexp_scale_exponent)
  call hdf5_read_attribute(file_handle,'PartType0/Temperature/h-scale-exponent'   ,Conv%Temperature%h_scale_exponent)
  call hdf5_read_attribute(file_handle,'PartType0/Mass/CGSConversionFactor'       ,Conv%Mass%CGSConversionFactor)
  call hdf5_read_attribute(file_handle,'PartType0/Mass/aexp-scale-exponent'       ,Conv%Mass%aexp_scale_exponent)
  call hdf5_read_attribute(file_handle,'PartType0/Mass/h-scale-exponent'          ,Conv%Mass%h_scale_exponent)
  call hdf5_read_attribute(file_handle,'PartType0/Density/CGSConversionFactor'    ,Conv%Density%CGSConversionFactor)  
  call hdf5_read_attribute(file_handle,'PartType0/Density/aexp-scale-exponent'    ,Conv%Density%aexp_scale_exponent)
  call hdf5_read_attribute(file_handle,'PartType0/Density/h-scale-exponent'       ,Conv%Density%h_scale_exponent)
  call hdf5_close_file(file_handle)

end subroutine read_cgs_conversions
